Pairing renormalization and regularization within the local density approximation 
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We discuss methods used in mean-field theories to treat pairing correlations within the local 
density approximation. Pairing renormalization and regularization procedures are compared in 
spherical and deformed nuclei. Both prescriptions give fairly similar results, although the theoretical 
motivation, simplicity, and stability of the regularization procedure makes it a method of choice for 
future applications. 
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I. INTRODUCTION 

One of the main goals of the low energy nuclear the- 
ory is to build a comprehensive microscopic framework, 
in which nuclear bulk properties, excitations, and low 
energy reactions can be described. For medium-mass 
and heavy nuclei, self-consistent methods based on the 
Density Functional Theory (DFT) 0, 0] have already 
achieved a level of precision that allows for analysis of 
experimental data for a wide range of properties of nu- 
clei throughout the chart of the nuclides. For example, 
the self-consistent Hartrcc-Fock-Bogoliubov (HFB) mod- 
els based on the Skyrme energy functionals 0, 0, Q are 
able nowadays to reproduce nuclear masses with an rms 
error of about 700 keV [|| 0] • The development of a uni- 
versal nuclear density functional, however, still requires 
a better understanding and improved description of the 
density dependence, isospin effects, pairing force, many- 
body correlations, and symmetry restoration. 

Nuclear pairing is an important ingredient of the nu- 
clear density functional, and it becomes crucial for open 
shell nuclei, in particular weakly bound systems, where 
the effects of coupling to continuum become significant 
In this case, the BCS model is not adequate Q 
and the fully self-consistent HFB approach must be used. 

In most HFB applications, pairing interaction is as- 
sumed to be either in the form of the finite range Gogny 
force or the zero-range, possibly density-dependent, 
delta force @, H E] • Gogny interaction in the pairing 
channel can be viewed as a regularized contact interac- 
tion, with regularization fixed through the finite range. 
The resulting pairing field is, however, nonlocal. 

Calculations using the contact interaction are numer- 
ically simpler, but one has to apply a cutoff procedure 
within a given space of single-particle (s.p.) states pHT^. 
When the dimension of this space increases, the pairing 
gap diverges for any given strength of the interaction. 
Therefore, the pairing strength has to be readjusted for 



each s.p. space. Thus the energy cutoff and the pairing 
strength together define the pairing interaction, and this 
definition can be understood as a phenomenological in- 
troduction of finite range d Such a procedure is 
usually referred to as the renormalization of the contact 
pairing force. It is performed in the spirit of the effective 
field theory, whereupon contact interactions are used to 
describe low energy phenomena while the coupling con- 
stants are readjusted for any given energy cutoff to take 
into account neglected high energy effects. 

The renormalization procedure for the zero-range pair- 
ing interactions has been explored in Ref. [8j using the 
numerical solutions of the HFB equations. It has been 
shown that by renormalizing the pairing strength for each 
value of the cutoff energy one practically eliminates the 
dependence of the HFB energy on the cutoff parameter. 

Recently, the issue of contact pairin g fo rce has been 
addressed in Refs. H E E lH EM E E E3, 
suggesting that the renormalization procedure can be re- 
placed by a regularization scheme which removes the cut- 
off energy dependence of the pairing strength. In subse- 
quent papers, this regularization scheme has been applied 
to properties of the infinite nuclear matter 119 1 . sphe rical 
nuclei 0, E] , and trapped fermionic atoms |loll24| . 

In this study we investigate the stability of the regular- 
ization scheme with respect to the cutoff energy for both 
spherical and deformed nuclei. Differences between the 
HFB results emerging from the pairing renormalization 
and pairing regularization procedures are analyzed. 

The HFB and Skyrme HFB formalisms have been ex- 
plained in great detail in many papers (see, e.g., Refs. 
|l2ll25j ]'). The notation used in t he p resent paper is con- 
sistent with that of Refs. 0, l2»| . This work is orga- 
nized as follows. Sec. Ull gives a brief introduction to the 
pairing renormalization and regularization schemes. In 
Sect. ITTT1 we explain the numerical framework used. The 
comparison between pairing regularization and renormal- 
ization techniques, studied for a large set of spherical and 
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deformed nuclei, is discussed in Sec. II VI Finally, the sum- 
mary and conclusions are given in Sec. [3 



II. THE CUTOFF PROCEDURES 

A. Pairing Renormalization Procedure 

Within the HFB theory, the energy cutoff can be ap- 
plied either to the s.p. or to the quasiparticle spectrum. 
The first option is used when the HFB equations are 
solved within a restricted s.p. space. However, the s.p. 
energies play only an auxiliary role in the HFB method, 
and the cutoff applied to the quasiparticle spectrum is 
more justified. This is done by using the so-called equiv- 
alent s.p. spectrum |l2^ : 



e„ = (1 - 2P n )E n + p, 



(2.1) 



where E n is the quasiparticle energy and P n denotes the 
norm of the lower component of the HFB wave function. 

Due to the similarity between e n and the s.p. energies, 
one takes into account only those quasiparticle states for 
which e n is less than the assumed cutoff energy e cut . 

It was shown || that for a fixed pairing strength the 
pairing energies depend significantly on the energy cut- 
off. Within the renormalization scheme employed in this 
work, we use the prescription of adjusting the pairing 
strength to obtain a fixed average neutron pairing gap 

A = - — / d 3 rdV H™, rV)p(rV, ra), (2.2) 

in 120 Sn equal to the experimental value of 1.245 MeV. 
In Eq. H2.2J) N is the number of particles, p is the par- 
ticle density, and h is the pairing Hamiltonian (see Ap- 
pendix . 

Such a procedure almost eliminates the dependence of 
the HFB energy on the cutoff 



B. Pairing Regularization Procedure 

Using the HFB equations and properties of the Bo- 
goliubov transformation (see appendixlAlfor details), one 
concludes that the local abnormal density p has a singular 
behavior when e cu t — ► oo. The standard regularization 
technique is to remove the divergent part and define the 
regularized local abnormal density /5 r (r) as 

p r (r) = lim [p(r - x/2, r + x/2) - /(r, x)] , (2.3) 

x— »0 



where / is a regulator which removes the divergence at 
x = 0. 

For cutoff energies high enough, one can explicitly 
identify 0, 0, components generating diver gen ce 
in the abnormal density (see, e.g., Eq. (21) of Ref. 



/(r,x) = «^M£) + i|> G ,.( r + x/,, r -x/ 2 ), 

(2.4) 

where G M is the s.p. Green's function at the Fermi level 
p in the truncated space, M* is the effective mass, and 
the Fermi momentum is 



Mr) = ^r^V^Wl (2-5) 



with U being the self-consistent mean-field potential. 

The first term in Eq. (|2.4|l comes from the MacLau- 
rin expansion with respect to x; it guarantees that the 
regularization procedure does not introduce any constant 
term to the abnormal density and that /(r, x) solely rep- 
resents the divergent part of p. 

Using the Thomas-Fermi approximation, the local s.p. 
Green's function G M (r) := Gpfr, r) becomes [TsL ITflj 



I i-fccut(r) 

G u (r) = ——, r lim / , ,, , 



k 2 dk 



^ 2 ^Jo p-^-U(r) +i j 
where the cutoff momentum is given by 
\/2M*(r) 



(2.6) 



fccut(r) 



\Acut+^-E/(r). (2.7) 



The regularized pairing Hamiltonian and the pairing 
energy density may be written, respectively, as [lSj : 

h(r) = g(r)p r (r) = g eff (r)p(r) (2.8) 



Wp«>(r) = -5e//(r)p(r) 2 , 



(2.9) 



where the effective pairing strength 0, 0, ^ , 



9W 2 



after calculating integral l|2.6|) . can be expressed in the 
form: 
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9eff(r) 



_L M*(r)fc cut (r) 

g(r) 2tt 2 H 2 



( 1 _ kp ( : 
V 2fe cut 



M*(r)fc cut (r) 

27T 2 R 2 



r ) l n fcout(r)+fc F (r) 
(r) 111 feout(r)-fcjr(r) 



iMl)l arctan ^4)l 

fecut(r) fccut(r) 



k F (r) 2 > 



fc F (r) 2 < 



(2.11) 



In this regularization scheme, only the Green's func- 
tion is calculated using the Thomas-Fermi approxima- 
tion. The densities, potentials, and chemical potential 
are determined self-consistently within the HFB theory. 
Consequently, the Fermi momentum (|2.5(l depends on mi- 
croscopic HFB quantities. According to the sign of k F , 
one of the exp ressions i|2.11|) is used. 

In Ref. [23 a different regularization scheme has been 
proposed that involves truncation below and above the 
Fermi level. However, the HFB calculations in the quasi- 
particle basis should be performed for a high cutoff en- 
ergy of 50MeV and higher ||- Since the magnitude of 
the self-consistent mean field U is also about 50 MeV, for 
such a high cutoff energy both methods are equivalent. 
The Thomas- Fermi approximation requires that, in order 
to obtain results independent of e cut , its value should be 
high enough for k cut to be real everywhere. 

Through the density dependence of g e ff, k cu t, and k F , 
there appear rearrangement terms in the self-consistent 
mean-field potential: 

S7~Cpair &g e f f n ~2 

/ = ~T^P =P x 
Sp Sp 



x ( 9 9eff Sg 
V dg Sp 



dg eff Sk F 
dkp Sp 



dgeff Sk cut * 
dk cut Sp , 



.(2.12) 



The first term in Eq. I|2.12|l is similar to the usual re- 
arrangement term, while the other two terms associated 
with the regularization procedure are entirely new. It is 
easy to check that all the terms appearing in Eq. (|2.12|) 
are continuous at the classical turning point k F (r) = 0. 

In Eq. H2.9|) . the pairing energy density is divergent 
with respect to the cutoff energy. However, the pairing 
energy itself is not an observable, and in order for the 
energy density functional to be independent of the cut- 
off, other terms have to cancel out this divergence. As 
discussed in Refs. 0,0,0,0, the kinetic energy den- 
sity t has the same type of divergence as the abnormal 
density p, and the sum 



■in+pair 



(r) 



2M*(r) 



r+7;9eff(r)p 2 (r) (2.13) 



does converge. 

Various contributions to the total HFB energy as func- 
tions of the cutoff energy are shown in Fig. ^ The 
total energy is stable with respect to e cut , although 
some of the components of the total energy vary signifi- 
cantly. As expected from Eq. I|2.13|l . two terms exhibit- 
ing large fluctuations are the kinetic term (with vari- 
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FIG. 1: Various contributions to the HFB energy for 120 Sn as 
a function of e cut . Calculations are performed using the SLy4 
Skyrme functional and mixed pairing interaction ffi.ip . 



ations of about 2 MeV) and pairing term (with varia- 
tions of about 1.3 MeV). Also, the momentum-dependent 
spin-orbit term, Es.0.1 has significant variations of about 
1 MeV. On the other hand, Skyrme and Coulomb ener- 
gies are fairly stable with respect to e cu t. 



III. NUMERICAL IMPLEMENTATION 

A. Numerical Framework 

As the pairing renormalization and regularization pro- 
cedures remove the divergent part of the abnormal den- 
sity in a different way, one can expect some numerical 
differences between both methods. In order to compare 
their results, we have performed numerical calculations 
using two numerical codes solving the HFB equations: 

• HFBRAD [23 - solves the HFB equations in the 
spherically symmetric coordinate basis. The max- 
imum angular momenta used in calculations were 
Jrnax = 39/2 for neutrons and j ma x = 25/2 for pro- 
tons. 
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FIG. 2: Total energy (top) and neutron pairing gap (bottom) 
in 120 Sn without (left) and with (right) pairing renormaliza- 
tion applied. Results are shown for volume (gray) and mixed 
(black) pairing. The total energy is plotted relative to the 
values obtained for the cutoff energy of e C ut=60 MeV. 



• HFBTHO |26| - diagonalizes the HFB problem in 
the axially symmetric transformed harmonic oscil- 
lator (HO) basis. Unless stated otherwise, we use 
N osc — 20 HO shells in the basis. 

In our calculations, we use the SLy4 and SkP 0] 
parameterizations of the Skyrme functional in the p-h 
channel and the contact density-dependent force in the 
p-p channel, which leads to the pairing energy density of 
the form: 



Impair (*") 



1 



'(r)/5(r) 2 = 



1* 



1-Vi 



Po 



p(r) 2 , (3.1) 



where po = 0.16 fm~ 3 . For V\ = the resulting pair- 
ing interaction is called volume pairing, while V\ = 1/2 
corres pon ds to the so-called mixed pairing prescription 
(Ref. [29| and references quoted therein). 



B. Pairing Renormalization 

Figure|3illustrates the importance of the pairing renor- 
malization procedure in the case of 120 Sn. Due to the 
constraint (|2.2|l on the pairing strength, the neutron av- 
erage pairing gap stays by definition constant, while the 
resulting total energy changes with the cutoff energy by 
a few hundred keV. On the other hand, without pairing 
renormalization applied, the total energy and the average 



FIG. 3: Total energy (top) and neutron pairing gap (bottom) 
in 120 Sn for the two values of N osc (left) or using two box sizes 
(right). Calculations were performed using the mixed pairing 
interaction. 



neutron gap vary significantly with increasing dimension 
of the quasiparticle space. In this case, the total energy 
changes by several MeV. 



C. Pairing Regularization 

The total energy and the average neutron pairing gap 
in 120 Sn are shown in Fig. Rafter applying the pairing 
regularization procedure. The pairing strength Vo is kept 
constant; it reproduces the neutron pairing gap for 120 Sn 
at the cutoff energy of e cu t =60 MeV. 

In the left panels of Fig. [3J we show results obtained in 
the HO basis, while the results from the solution of the 
HFB equations in coordinate space are displayed in the 
right panels. One can correlate the coordinate-space and 
HO representations by introducing an 'effective box size' 
R ~ \j2N osc hlmuj Q . Using this formula, the basis of 20 
HO shells corresponds to a box radius of about 14.5 fm. 
Figure |21 demonstrates that the regularization procedure 
is stable with respect to the cutoff energy. Moreover, 
one obtains reasonable results already for fairly low cut- 
off energies of about 40 MeV. The variations in the to- 
tal energy in coordinate-space calculations do not exceed 
40keV, while they are about 150 keV in the HO expan- 
sion. The latter number does not decrease significantly 
with iV osc . 

The differences in applying the pairing regularization 
procedure in the coordinate-space and HO calculations 
can be explained by the different way the quasiparticle 
space is expanded in both approaches. The particle den- 
sity p is defined by the lower components of the quasi- 
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FIG. 4: Total energy of spherical 120 Sn (top) and deformed 
110 Zr (bottom) obtained with pairing regularization (black 
lines) for mixed pairing (left) and volume pairing (right). The 
results obtained without the rearrangement terms resulting 
from the variation of fc cut and h,F in Eq. I|2.12|l are also shown 
using (gray lines). 



particle wave functions, which are localized within the 
nuclear interior. On the other hand, the abnormal den- 
sity is defined by the products of the upper and lower 
components of the quasiparticle wave function. For the 
quasiparticle energies that are greater than the modu- 
lus of the chemical potential, the upper components of 
the quasiparticle wave function are not localized. There- 
fore, contrary to the normal density, the abnormal den- 
sity strongly depends on the completeness of the s.p. ba- 
sis outside the nuclear interior. 

In the coordinate-space calculations, the box boundary 
conditions provide discretization of the spectrum for the 
quasiparticle continuum states that are not localized. On 
the other hand, all the HO basis states are localized. Re- 
sults of stability with respect to the cutoff energy for the 
coordinate-space and HO calculations are, therefore, dif- 
ferent. As far as the description of nonlocalized states is 
concerned, the coordinate-space method is superior over 
the HO expansion method. 

Fluctuations in the total energy shown in Fig. [3] co- 
incide with 2j + 1-folded degenerate angular-momentum 
multiplets of states in spherical nuclei that enter the pair- 
ing window with increasing cutoff energy. This can be 
confirmed by performing a similar analysis for a deformed 
nucleus where the magnetic degeneracy is lifted. Such re- 
sults are shown in Fig.0]for deformed 110 Zr in comparison 
with spherical 120 Sn. One can see that the fluctuations 
of the total energy in 110 Zr are down to about 40keV. 

The steep increase of the total energy at the cutoff en- 
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FIG. 5: Two pairing regularization schemes applied to the 
case of 120 Sn: the Thomas- Fermi approximation [l9l ] (black 
line) and the free particle Green's function [l3| (gray line). 
Coordinate-space calculations were performed in a 15 fm box. 



ergies below 30 MeV results from neglecting quasiparticle 
states with significant occupation probability. This effect 
is more severe for the mixed-pairing than for volume- 
pairing calculations due to the surface-peaked character 
of mixed pairing fields. On the other hand, the stability 
with respect to the cutoff energy is similar in both cases. 

We have also tested the importance of the rearrange- 
ment terms arising as a result of the regularization pro- 
cedure. The gray lines in Fig. 0] show results obtained 
without taking into account the second and third term 
of Eq. H2.12fl . These terms lead to changes in the total 
energy of a few keV and can be safely neglected. 

Finally, we have tested the Thomas-Fermi approxima- 
tion used in the pairing regularization procedure. In- 
stead of adopting the Thomas-Fermi ansatz, one can per- 
form regularization using the free particle Green's func- 
tion [13J. As illustrated in Fig. [21 the convergence of the 
latter method is very slow; the Thomas-Fermi method is 
clearly superior. 



D. A link between the pairing renormalization and 
regularization procedures 

The renormalized and regularized pairing calculations 
are based, in fact, on two different effective interactions. 
Consequently, their results should be comparable only as 
much as their effective pairing strengths g e // are similar. 
By expanding Eq. (|2.11[) at very high cutoff energies 
(kp/k cu t << 1) 7 one obtains: 



9eff(r) 



1 



M*(r)g(r) 
2ir 2 h 2 



fccut (r) g{r), (3.2) 



which has the form of 3 e // = a 9- For the volume pairing, 
the proportionality factor a is p-dependent only through 
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FIG. 6: Ratio between the effective pairing strengths for pair- 
ing regularization and renormalization, g^f/g^ff, for the vol- 
ume (upper panel) and mixed pairing (lower panel) in 120 Sn 
for several values of e cut . 



the weak density dependence of the effective mass M*. 
On the other hand, for the mixed pairing, it also depends 
on p through the density dependence of g. Therefore, 
while for the volume pairing the renormalization proce- 
dure may be considered as a fair approximation to the 
regularization scheme, this is not the case for the mixed 
pairing, or - more generally - for any density-dependent 
pairing. Still, this approximate equality of the effective 
pairing strengths for the pairing regularization and renor- 
malization is an explanation of the remarkable stability 
of the total energy in phenomenological pairing renormal- 
ization treatment (see Fig. |5J) , and it also explains why 
results obtained for the volume pairing are more stable 
than those in the mixed pairing variant. 

This effect can be clearly seen in Fig. The ratio 
between the effective pairing strengths in the regulariza- 
tion and renormalization methods is much closer to unity 
for the volume pairing than for the mixed pairing in the 
region of space, where the pairing energy density is max- 
imal. 



IV. COMPARISON BETWEEN PAIRING 
RENORMALIZATION AND REGULARIZATION 
PROCEDURES 

In this section, we present a comparison between pair- 
ing renormalization and regularization procedures ap- 



plied to a large number of nuclei. As representative re- 
sults, we discuss results obtained for the drip-to-drip line 
isotopic chains of spherical Sn nuclei as well as for de- 
formed Dy nuclei. Calculations are performed for the 
volume and/or mixed pairing interactions by using the 
HFB+THO approach. 



A. Spherical Nuclei 

Figure[7]displays differences between the pairing renor- 
malization and regularization procedures for the Sn iso- 
topes. Calculations are performed with both volume and 
mixed pairing interactions. For the two-neutron sep- 
aration energies, the maximum difference between the 
renormalization and regularization schemes is about 100 
(300) keV for the volume (mixed) pairing. In the neutron 
gap, the corresponding difference is about 50 (100) keV, 
and in nuclear radii (not displayed) it is practically neg- 
ligible (about 0.01 fm). The largest differences show up 
in pairing energies - about 1 (3) MeV for the volume 
(mixed) pairing; however, total energy differences are 
much smaller - about 400 (800) keV. 

Analyzing the total energies obtained in both meth- 
ods, Fig. EJa), one can see that the pairing renormaliza- 
tion procedure gives systematically more binding. The 
differences are negligible for stable nuclei and nuclei near 
the proton drip line. They increase in mid-shell nuclei 
near the two-neutron drip line where the pairing effects 
are the largest, and then decrease towards the closed- 
shell nucleus 176 Sn located just at the two-neutron drip 
line. In general, both procedures give more similar re- 
sults in the case of volume pairing than in the case of 
mixed pairing. 

Recently, the pairing regularization procedure has been 
analyzed in the context of relativistic mean-held approx- 
imation [2^. In order to simulate the finite range con- 
tribution to the nuclear matter pairing gap coming from 
the Gogny pairing force, it was necessary to introduce 
strong density dependence in the pairing strength of the 
contact interaction. 

Using the regularization procedure and calculating the 
Sn chain with both volume and newly constructed (sur- 
face) contact interaction, the authors of Ref. [2^ have 
found differences in pairing energies of the order of 
20 MeV in the neutron-rich nuclei around 1 Sn. In our 
work, for the same nuclei, the differences in pairing ener- 
gies between volume and mixed pairing variants do not 
exceed 2.6 MeV. This comparison shows that the density- 
dependent contact interaction proposed in Ref. [23| is 
questionable for finite nuclei, despite its agreement with 
the finite-range Gogny pairing force in the infinite nuclear 
matter. 
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FIG. 7: Differences between pairing renormalization (RN) and regularization (RG) procedures for (a) total energies and 
neutron pairing energies, (b) two neutron separation energies, and (c) the average neutron gaps. The HFB+THO calculations 
are performed for the chain of the spherical Sn isotopes using the SkP Skyrme parameterization. 



B. Deformed Nuclei 

We applied the pairing renormalization and regulariza- 
tion procedures for the chain of deformed Dy isotopes. 
Differences between both sets of results are shown in 
Fig. El We show only the results with the mixed pairing, 
since, as in the spherical nuclei, the differences between 
both procedures are larger in this case. 

As seen in Fig. Etc), most of the nuclei considered 
are well deformed, and the deformations are practically 
the same within both procedures. Despite the fact that 
the maximum difference in the pairing energy is around 
3MeV (not shown), other quantities are very similar. 
The maximum difference in the total energy is about 
360 keV, in the two-neutron separation energy - 160 keV, 
in the pairing gaps - 110 keV, and in the rms radii (not 
shown) - less than 0.005 fm. 

V. SUMMARY AND CONCLUSIONS 

In this work, we investigated the pairing regularization 
method using the s.p. Green's function in the Thomas- 
Fermi approximation and found it to be very suitable 
for a description of spherical and deformed nuclei. We 
checked the stability of the method with respect to the 
cutoff energy and found fluctuations in the total energy 
below 200 keV. Fluctuations coming from the method it- 



self do not exceed 50 keV for the cutoff energy as low as 
30 MeV. However, if a still lower cutoff energy is assumed, 
the Thomas-Fermi approximation to the s.p. Green's 
function may no longer be valid. 

We found that the differences between pairing renor- 
malization and regularization procedures for volume and 
mixed pairing are rather small. Therefore, we conclude 
that physical conclusions previously obtained within the 
pairing renormalization scheme remain valid. Neverthe- 
less, we believe that the theoretical motivation and sim- 
plicity of the regularization method is preferred to a phe- 
nomenological renormalization scheme. 
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FIG. 8: Similar to Fig. Q except for the deformed Dy isotopes. Quadrupole deformations are displayed in panel (c). Mixed 
pairing interaction was used. 



APPENDIX A: DIVERGENCE IN THE 
ABNORMAL DENSITY 

In the DFT-HFB approach, the starting point is the 
Energy Density Functional (EDF) H[p, p], where p is the 
particle density and p is the abnormal density: 



fflTl"r2IT2 T 2 



I*) 



p(r 2 CT 2 T2,ricriTi) = ($1^ 
p(r 2 (J 2 T 2 ,Y 1 aiTi) = -2cri($|a ri _ (Tiri a r 



(Al) 
|$),(A2) 



where a and are the particle annihilation and creation 
operators, respectively, and |<I>) is the HFB state. In the 
following, we assume that |<I>) is a product of the neutron 
and proton states, l^?,,)!^). Therefore, the neutron and 
proton wave functions are not coupled, and in the nota- 
tion below we can, for simplicity, omit the isospin index 
with the understanding that all equations are separately 
valid for neutrons and protons. 

For the HFB state |$), the particle and abnormal den- 
sities can be written as [l2| : 

p(r 2 CT 2 ,ricri) = ¥>2i(r 2 CT2)^ 2 'i( r i cr i)) ( A 3) 

Ei>0 

p(r 2 cr 2 ,ricri) = - ^ <P2i(r2<T2)<Pu(Ti<7i)> ( A4 ) 

Ei>0 

where the two-component quasiparticle wave function (p 



is the solution of the HFB equation: 



M^^ri^i) h(r 2 a 2l r 1 (7 1 ) 



¥>ii(ricri) 
<£>2i(ri0i) 



h(r 2 a 2l r 1 ai 

= E, 



-/iu(r 2 cr 2 ,ri(Ti) 



Vii(r2CT2) 
^2j(r 2 cr2) 



(A5) 



for a given quasiparticle energy Ei. 

The HFB equations are a result of variational mini- 
mization of the energy density functional TC[p, p] with the 
constraint of the mean value of particles kept constant: 



Sn \(N) = 



{N)=N 



0. 



(A6) 



This condition defines the s.p. Hamiltonian and the 
pairing Hamiltonian h in the HFB equations l|A5(l : 



h tl (v 2 <7 2 ,r 1 (j 1 

= -V 



SH[p,p] 



<5p(ri<7i,r 2 c72 

h 2 



ra 2M*(r 2 <7 2 ,riCTi) 

5H[p,p] 



/i(r 2 cr 2 ,ricri) = 



(A7) 

V ri + U(r 2 a 2 ,r 1 ai) - p 
(A8) 



<5/5(ri£r 2 , r 2 a 2 ) 



where M* is the effective mass and U is the self-consistent 
mean-field potential. In the following derivations, the 
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spin-orbit term is omitted as unimportant in the regular- 
ization scheme, although it is, of course, always included 
in calculations. 



By multiplying the HFB equations (|A5|) by vector 
[v?2i> — ViiL integrating over coordinates and summing 
over all the positive energy HFB solutions, one obtains: 



53 EiJ d 3 r 2 [ <^(r 2 cr 2 ), -<^(r 2 cr 2 ) ] 



Ei>0,a 2 



</>li(r 2 CT 2 ) 



53 // d 3 ri d 3 r 2 [ <p* 2l {v 2 a 2 ), -<^(r 2 a 2 ) ] 

E i >0,cr 1 a 2 ^ 



Mr2°" 2 ,rio-i) h(r 2 a 2 ,r 1 a 1 ) 
h(r 2 a 2 ,r 1 a 1 ) -/i M (r 2 cr 2 , rid) 



¥>2i(ri<ri) 



(A9) 



53 / rf 3 n {V2i( r lo-l)Vli( r io-l) - ^i(ricri)y>2i(ricri)} = 

d 3 v x d 3 v 2 \ ipl i (r 2 a 2 )h fl (r 2 a 2 ,r 1 a 1 )ipi i (r 1 ai) + <^ 2l {r 2 a 2 )h(Y 2 a 2 , v\(J\)tpn(T\(T\) 



Ei>0,cr 



E 



+<pl i (r 2 a 2 )h ll (r 2 <j 2 ,r 1 <7i)(p 2i (r 1 ai) - ^* u {r 2 <7 2 )h(v 2 a 2 ,v 1 a 1 )ipi l {r 1 ai) 



(A10) 



Since for every HFB solution ([ipn, (p 2 i\, Ei) there exists also an orthogonal solution ([ip2i, —'fin], —Ei), the left-hand 
side of Eq. (|A10|) vanishes as a sum over scalar products of orthogonal wave functions. 
For local and spin-independent Hamiltonians and h, Eqs. (|A8|) and (| A8|> read 



^(r2C 2 ,ricri) = -V r 



ft(r 2 <T2,ri(Ti) = h(r 2 )S(r 2 - r 1 )S a2 ^ 1 



8(r 2 - ri)<5 (T2)CTl Vr 1 + (U(r 2 ) - n)6(r 2 - r 1 )S l72t(Tl , 



(All) 
(A12) 



Note that for an attractive pairing force, the local pairing potential h(r) = — A(r) is negative, where A(r) is the 
standard position-dependent pairing gap. By defining function J- £cut as 



53 [^«( r 2^Mi(>iO-) +^2i(r2cr)^2i ( r ic)] = ^ecutfo ~ ri). 



(A13) 



Ei>0,CT 

and using expression l|A4(l for the abnormal density, one obtains after integrating the kinetic-energy term by parts: 

h 2 







d i Y 1 d i v 2 5(r 2 - r x ) 



- I d 3 rd 3 x<5(x) 



h(r 2 )[F €cM (r 2 -r 1 )-2p{r 2 ,r 1 )} 

h 2 (1 



2M 



-V r2 V ri + C/(r 2 )-/i 2 / 5(r 1 ,r 2 



h(r) [^ eeut (x)-2p(r,r)] + 



2M* \4 



- Vi + C/(r) - M 2p(r - x/2, r + x/2) 



(A14) 



r 



where 



r = 



ri + r 2 



When the summation over positive quasiparticle energies 
is extended to infinity, the completeness relation implies 
(A15) that 



and 



x = r 2 ri, 

/o(r 2 ,ri) = 53/ 3 ( r 2< J ,ricr), 
/5(r 2 ,ri) = 53^( r 2 cr > r i cr )- 



(A16) 



^e cut (r 2 -T X ) =S(r 2 - ri), 



(A19) 



(A17) and the only term in Eq. I|A14|I capable of canceling out 
this singularity is Vx/5(r — x/2, r + x/2). Therefore, the 

(A18) Laplacian of the abnormal density Vx/5(r — x/2, r + x/2) 
must be singular at x = 0. Moreover, using the expres- 
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= -AnS(r), 



(A20) 



it is clear that due to the zero-range pairing interaction 
abnormal density p has an ultraviolet l/x divergence: 



p(r - x/2, r + x/2) 



h{v)M*{v) 
4tt?i 2 |x| 



(A21) 



x^O 
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